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1 .  Introduction 

We  consider  iterative  methods  for  solving  the  linear  comple¬ 
mentarity  problem  (LCP) : 

Find  X  in  :  x^O,  Mx+q^O,  x^(Mx+q)  =0  (1) 

where  M  is  a  given  n-by-n  real  matrix  and  q  is  a  givcu  real  n- 
vector.  The  general  LCP  is  NP-complete  (Ref.  1)  and  therefore  very 
difficult  to  solve  even  for  problems  with  moderate  size,  say  n  =  50. 
However,  in  many  applications  the  matrix  M  has  some  nice  proper¬ 
ties,  e.g.  M  is  positive  semidefinite  or  all  its  principal  minors 
are  positive  (i.e.  M  is  a  P-matrix) ,  and  the  problem  can  be  solved 
by  some  direct  pivoting  methods,  e.g.  the  principal  pivoting 
method  (Ref.  2)  or  Lemke's  method  (Ref.  3),  or  a  semi-iterative 
pivoting  method  by  Shiau  (Ref.  4). 

Although  a  pivoting  method  gives  a  very  accurate  solution  in 
a  finite  number  of  steps  ,  i*:  ir  ^.y  not  be  suitable  for  some  large 
sparse  problems.  Since  pivoting  can  easily  destroy  the  sparsity 
after  a  number  of  steps,  0(n  )  storage  space  is  needed  for  the 
pivoting  even  if  M  has  only  0(n)  nonzero  entries.  Iterative 
methods  that  preserve  sparsity  and  work  only  on  the  nonzero  entries 
of  the  input  data  are  much  more  attractive  for  very  large  prob¬ 
lems,  for  the  storage  requirement  is  greatly  reduced  and  the 
processing  could  be  faster. 
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Many  iterative  methods  have  been  proposed  (e.g.  Refs.  5-9). 
However,  they  were  developed  primarily  for  the  case  that  M  has 
some  kind  of  strict  positivities .  For  example,  Mangasarian ' s 
algorithm  (Ref.  9)  converges  when  M  is  symmetric  and  strictly  co¬ 
positive,  or  copositive-plus  and  there  exists  an  x  such  that 
MX+q>0.  In  one  of  the  most  important  applications  where  the  LCP 
is  formulated  for  solving  the  underlying  quadratic  program  (Ref.  2), 
the  matrix  M  is  only  positive  semidefinite  and  not  definite.  (There 
are  special  cases,  e.g.  the  quadratic  program  is  strictly  convex 
and  separable,  that  existing  methods  may  apply  (e.g.  Refs.  10-12)). 
In  this  work,  we  present  a  sparsity-preserving  iterative  method 
that  can  solve  nons5rmmetric  LCP's  when  M  is  only  positive  semi¬ 
definite  . 

We  give  a  brief  word  about  notation.  The  n-dimensional 
Euclidean  space  is  denoted  by  R"^  .  Vectors  in  R”  are  column  vectors 
and  denoted  by  x,  y,  z,  etc.  Superscripts  are  used  to  denote 
different  vectors  such  as  x° ,  ,  but  the  superscript  T  denotes 

the  transpose,  e.g.  x  ,  (x  )  ,  M  . 

2 .  Method  and  Convergence  Theorem 

Let  us  consider  first  the  following  quadratic  program: 


0  =  min  f(x)  subject  to  xeS , 


(2) 
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where  f  is  a  convex  quadratic  function,  S  is  defined  by  a  finite 
number  of  linear  inequalities,  and  0  =  f(x)^f(x)  for  some  x  in 
S,  for  all  X  in  S .  The  intension  is  to  set  S={xix^0,  Mx+q^O) 
and  f(x)  =  x’^(Mx+q)  ,  then  all  the  assumptions  of  f  and  S  are 
satisfied  provided  that  M  is  positive  semidefinite  and  the  LCP 
has  a  solution.  In  this  setting,  every  solution  of  (2)  is  a 
solution  of  (1)  and  vice  versa.  We  present  the  following  iter¬ 
ative  scheme  and  its  convergence  theorem.  The  proof  of  a  more 
general  scheme  can  be  found  in  (Ref.  4)  and  is  very  similar  to 
that  In  (Ref.  13).  We  include  the  proof  here  to  make  the  paper 
self-contained. 


2.1  Iterative  Scheme 


Step  1. 


Ic  Ic 

Given  x  e  S,  stop  if  f(x  )  =  0  (within  a  fixed  tolerance). 
Otherwise,  find  y  satisfying 


y  e  S,  f(x'^)  +  vf(x’")^(y’^  -  x^)<0 


(3) 


Step  2. 


Let  p'^  =  y*^  ~  be  the  direction  of  search, 

j  k+l  k  ,  _ 

find  X  =  X  +  t,  p 

k 

where  t  is  defined  by 

iC 


K  ]r 

t  =  arg  minimize  f(x  +  tp  ) 
^  0<  t<  1 


(4) 


Step  3.  Go  to  Step  1. 
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The  computation  of  y  is  discussed  in  the  next  section.  Since  f(x) 
is  quadratic,  the  objective  function  of  (4)  is  a  quadratic  function 
of  t,  in  fact, 

fCx’"  +  tp^)  =  f(x*")  +  (vf(x’")^p^)t  +  %  t2(p’")'^7^f(x^)p^ 

=  :  a^^  +  bj^t  +  %  Cj^t^  (5) 

Note  that  b^  -f(x  )<0  by  (3),  and  that  c,^^  0  by  the  convexity  of 
f(x).  It  is  very  easy  to  compute  t^^  as  shown  by  the  following 
lemma . 


Lemma  2.1.  Let  g(t)  =  a  +  bt  +  %  ct"^  ,  where  b<o,  c20.  Then 
t  ;=  /-b/c,  if  c  > -b , 


1,  otherwise,  is  the  minimizer  of  the  following  problem, 
minimize  g(t),  (6a) 

O^t^l  (6b) 

moreover , 

g(0)  -  g(t) >  -  ^b  t .  (7) 


Proof  the  case  c  =  o  is  trival.  Suppose  c>0. 


So  g'(t)<0  for  0<t<-b/c,  and  g'(t)>0  for  t>-b/c.  Therefore  g(t) 

is  strictly  decreasing  on  [0,-b/c]  and  increasing  on  [-b/c,+  «=]  .  If 
-b/c  <1,  i.e.  c^-b,  then  t  =  -b/c  solves  (6)  and  it  is  straightforward 

to  see  that  (7)  holds  as  an  equality.  If  -b/c  >1,  then  t  =  1  is  the 

constrained  minimum  and 

g(0)-g(l)  =  -b-%  c> -b/2 


proving  ( 7 ) . 
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Theorem  2.1.  Convergence  Theorem.  Assume  that  f(x)  is  convex  and 

quadratic,  S  is  nonempty  and  0  =  min  f(x).  Let  {x'^,y’^}  be  generated 

xes 

•  k  k 

as  in  2 . 1 .  If  {y  }  is  bounded,  then  converges  monotonously 

to  zero. 

Proof  By  Lemma  2.1, 

f(x'^)  -  f(x’^'^^)>-  %  t^fCx’^)  (8) 

where  b^^  is  defined  in  (5),  and  the  last  inequality  follows  by  (3). 
Summing  up  (8)  for  k  =  0,  1,...,N,  we  have 

f(x°)  -  f(x^‘^^)>%  I  t  fCx’^) 

k=o  ^ 

Since  the  left-hand  side  is  bounded  for  all  N,  the  positive  series 
00 

^  t  f(x^)  converges.  Hence 

k=o 

lim  t  f(x^)  =  0  (9) 

k->-“  k 

By  Lemma  2.1,  is  either  1  or  tj^  =  By  assumption  {y^} 

k,  k  0 

is  bounded,  so  is  {x  }  since  x  is  a  convex  combination  of  x  and 

y  j  =  0,  l,...,k-l.  Therefore,  is  also  bounded. 
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So 

t  ^  min  {1,  b^/V}-min  {1,  f(x’")/V}  (10) 

where  V  is  a  bound  of  Since  f£(x^)}  is  monotonously 

decreasing  by  (8),  if  it  did  not  converge  to  zero,  then 
f(x^)  >  5  >  0  for  all  k,  for  some  <S .  Then,  by  (10), 
tj^f(x’^)  ^  min  {1,  f  (x’')/V}*f(x’")  >  min  {5,5^/V}, 

which  contradicts  to  (9).  ~  □ 

Corollary  2.1.  If  we  set  f(x)  =  x'^(Mx+q),  S  =  {xjxsO,  Mx+q>0}, 

where  M  is  positive  semidefinite .  Assume  S  is  nonempty  and  {y  } 

k  “ 

is  bounded,  then  lim  d(x  ,S)  =  0 

k^oo 

where  S  is  the  solution  set  of  (1),  and  d(x,S)  is  the  distance 

between  x  and  S  defined  by 

d(x,S)  =  inf  1 1  x-yjl 
y  eS 

Proof  It  follows  by  Theorem  2.1  and  Corollary  2.1  of  (Ref.  14), 
which  shews  that  for  xeS 

d(x,S)  .=  0(f  (X)  +  /f(x))  .  □ 

The  inequality  in  (3)  can  be  considered  as  a  cutting  plane 
which  cuts  off  part  of  S  including  the  current  point  x  ,  and  the 
remaining  part  still  contains  the  whole  S,  because  of  the 
convexity  of  f(x).  In  the  case  that  M  is  not  positive  semidefinite, 
but  has  other  positivity  properties  such  as  that  all  principal 
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minors  of  which  are  positive,  the  scheme  can  be  modified  by  a 
weaker  cutting  plane  which  cuts  off  less  than  the  one  in  (3) 
does,  so  as  to  guarantee  that  y*^  exists,  The  new  cut  is  defined 
by 

9j^f(x^)  +  Vf  (x’^)'^(x-x^)  <  0 

where  9j^  <  1.  The  number  0j^  should  be  as  large  as  is  allowed 
(for  the  existence  of  y  )  to  speed  up  the  algorithm.  For  sim¬ 
plicity  of  the  paper,  we  omit  this  case  and  concentrate  on  that 
M  is  positive  semidef  inite ,  thereby  9^^  is  always  chosen  to  be  1. 

Jr 

3 .  Solving  y  By  SOR 

The  main  effort  in  the  iterative  scheme  of  2.1  for  solving 

(I)  is  the  computation  of  y  satisfying  (3) .  This  is  to  find 
a  feasible  point,  in  the  k-th  iteration,  of  the  system  of 
linear  constraints: 

My  +  q  >  0,  y  >  0,  Vf  (x^)  (y-x^)  <  -f(x^)  (11) 

k  k- 1 

Note  that  X  and  y  satisfy  all  but  the  last  inequality,  i.e. 
the  cut,  and  therefore  they  can  be  selected  as  starting  points 
for  solving  y  .  There  are  a  few  non-pivoting  methods  for  solving 

(II)  which  involve  the  inversion  of  a  matrix  in  each  iteration. 

For  example,  by  replacing  the  last  inequality  by  equality,  one 
can  easily  employ  the  basic  (single  phase,  no  sliding  variable) 
Karmarkar  method  (Ref.  15)  to: 

Vf(x  )  X  -f(x  )  =  min  Vf(x  )  y 

s.t.  My  +  q  >  0,  y>0 
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We  propose  to  use  a  successive  overrelaxation  method  (SOR)  which 
needs  no  .atrix  inversion  and  is  sparsity-preserving  (Refs.  16-18). 
The  name  SOR  is  more  commonly  referred  to  the  relaxation  method 
for  solving  systems  of  n  linear  equations  in  n  unknowns  (e.g.  see 
Refs.  19-20).  The  SOR  for  system  (11)  of  linear  inequalities  is 
different  but  similar. 

3 . 1  SOR  for  Systems  of  Linear  Inequalities 

Problem  3.1;  Find  y  in  R^  satisfying  A^Y  ^  b,  i  =  1,  2,.,.,m, 
where  R"  ,  b^eR  are  given. 

Algorithm  3.1  Step  1.  Given  ,  j^O,  stop  if  z^  satisfies  (within 
a  tolerance)  all  the  inequalities.  Otherwise,  pick  the  next 
inequality  which  z^  violates,  say  A^Y<  b^ ,  compute 

z^  +  u((b.  -A'fz^)/A^A  )A,  (12) 

11  111 

where  0<u  <  2 . 

Step  2.  Go  to  Step  1. 

In  Step  1,  the  selection  of  the  next  violated  inequality  is 

based  on  the  cyclic  order  1,  2 . m,  1,  2,..  .  One  alternative 

is  to  choose  the  one  that  z  violates  most,  i.e.  choose  i  so 
that  (b^-A^z^)/l|  A^  j|  is  most  negative.  The  strict  cyclic  order 
can  be  relaxed  to  almost  cyclic  order  defined  as  follows.  Let 
i^  ,  i^  ,  .  .  .  ,  i  ,  ...  be  the  indices  of  the  inequalities  to  be  checked 
for  possible  violation  and  relaxation  (12),  the  order  is  called 
almost  cyclic  if  1  '  i .  ■ 


m ,  and 
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there  exists  an  integer  C  such  that  {1 ,2  , . . ,  ,m}E{i  ,i  ,..,i  } 

V  X/  ^  1  ^  C 

for  all  £^1.  The  strict  cyclic  order  is  almost  cyclic  with  C  =  m 
(Ref.  18).  When  m  =  1,  (12)  moves  along  -A^  to  to  satisfy  the 

i-th  constraint.  The  constraint  is  said  to  be  relaxed  in  the  process. 
It  is  an  overrelaxation  if  u  >1,  and  an  underrelaxation  if  u<l. 

Theorem  3.1 

Let  y  be  any  solution  of  the  Problem  3.1,  let  {z^}  be 
generated  by  the  SOR  method  (with  an  almost  cyclic  order  or  most- 
violation  rule). 

(a)  If  0<u<2,  then  lz^''’^-yj  <  jz^-yj,  and  {z^}  converges  to  a 
solution. 

(b)  If  there  exists  y  such  that  A.  Y <  b.  for  i  =  l,2,...,m, 

then  by  choosing  u  =  2,  the  algorithm  will  find  a  solution  in 
a  finite  number  of  iterations. 

Proof  See  (Refs.  16-18).  □ 

One  can  vary  u  in  each  iteration,  say,  choose  for  the  j-th 
iteration.  The  conclusion  of  (a)  remains  true  if  0<e^n^<2  -  e 
for  all  j  for  some  e. 

4 .  Computational  Details  and  Heuristic  Enhancements 

4 . 1  Since  SOR  is  an  infinite  procedure,  we  may  not  solve  y  exactly. 
However,  we  can  get  a  point  feasible  to  within  a  small  tolerance 
in  a  finite  number  of  iterations.  Since  any  solutions 
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of  the  LCP  (1)  satisfies  (11),  each  SOR  iteration  moves  the  point 
closer  to  the  solution  set  Theorem  3.1(a).  The  overall  solu¬ 
tion  process  is,  therefore,  applying  the  SOR  iteration  (12)  until 
the  movement  is  small,  followed  by  the  line-search  to  find  x 
and  updating  the  cutting  plane,  then  resuming  the  SOR  iteration. 

If  we  select  y’^  instead  of  as  the  starting  point  when  SOR 

restarts,  then  we  have  the  monotonously  decreasing  of  {d(z^,  S)} 
as  well  as  {f(x’^)},  where  {z^}  refers  to  all  the  points  generated 
by  SOR  which  contains  {y  }  as  a -subsequence . 

A . 2  When  applying  SOR,  the  n  nonnegativity  constraints  xs  0  can  be 
relaxed  easily  by  replacing  the  negative  components  of  z^  with 
zeroes  or  small  positive  niambers ,  for  overrelaxation.  Hence,  they 
can  be  relaxed  as  often  as  possible,  i.e.  every  time  any  component 
of  z^  becomes  negative  after  relaxing  one  of  the  rest  n+1  con¬ 
straints,  that  component  can  be  made  nonnegative.  The  convergence 
of  {z^}  still  holds,  for  the  resulting  order  is  almost  cyclic 
(defined  in  3.1  with  C  =  n(n+l)). 

A . 3  Heuristically ,  if  a  constraint  x^ ^  0  is  violated  very  often, 
we  can  set  x^=0,  thereby  reduce  the  problem  size.  The  reasoning 
is  that  if  {z^}  converges  to  x  and  x^  >  0  were  nonzero,  then 
should  also  be  positive  for  all  but  finitely  many  j.  Similarly,  if 
w.  >  0,  where  (w^ jW^ , • • • )^ :=  Mx+q ,  is  violated  very  often,  we  can  set 
"  =0.  i.e.  solve  (Mx+q)  =  0  for  some  x„  and  thereby  eliminate  var- 

i  1  X, 

iable  x  and  constraints  x  ^0,  w.^O.  This  clearly  requires  some 

1 

computation,  but  it  can  easily  be  justified  if  M  is  very  sparse. 

Since  we  also  relax  the  nonnegativity  constraints  of  x^^,  x^^  should 
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be  chosen  among  only  variables  for  which  the  current  values  are 
sufficiently  large.  It  is  possible  that  we  may  erroneously  set 
or  w^  to  zero,  but  we  can  easily  back  up  by  restoring  the  variables 
and  constraints.  For  example,  if  the  decreasing  of  f(x  )  becomes 
unsatisfactorily  slow  after  setting  w.  =0  and  removing  x  as  a 

X  X> 

consequence,  we  can  compute  the  value  of  x  by  the  equation  w,  =0, 
restore  the  constraints  w^^:  0,  x^>  0,  and  resume  the  computation  of 
Xj^  in  the  SOR  iteration. 

By  similar  argiiments  if  a  constraint  x^s  0(or  w^>0)  has  not  been 
violated  lately,  we  may  set  w^  =  0  (or  x^  =  0). 

Ic  + 1 

4 . 4  The  line -search  computes  x  as  the  minimizer  of  f(x)  on 

k  k 

the  line  segment  with  endpoints  x  and  y  .  We  may  do  higher  dimen¬ 
sional  searches,  e.g.  by  minimizing  f(x)  on  the  triangle  with  ver¬ 
tices  x^ ,  y^~^,  and  y^.  This  is  not  much  harder  than  the  line" 
search  since  f(x)  is  quadratic. 

4 . 5  Strictly  speaking,  y  does  not  satisfy  (11)  exactly,  rather, 
it  satisfies: 

y’^eSj^  :=  {y|My  +  q  +  e|^e>0,  y  +  Gj^e>0}  ,vf  (x’^)'^(y’^-  x’^)<  -f  (x’^) 

where  e  >  0  is  the  tolerance  in  iteration  k,  and  e  is  the  vector 
k 

of  ones.  However,  after  replacing  S  in  (3)  by  ,  Convergence 

Theorem  2.1  remains  valid,  i.e.,  we  still  have  f(x’^);0.  Corollary 

2.1  needs  to  be  modified.  If  we  choose  so  that  Gj^-*-0,we  can  still 

have  lim  d(x'^,S)  =  0 ,  by  Theorem  2.7  of  (Ref.  14). 
k-*-“ 
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5.  Numerical  Results 

A  computer  program  is  written  on  an  IBM  personal  computer 
(PC  AT)  for  a  special  class  of  LCPs  for  which  the  matrix  is  of 
the  form: 


where  is  either  1  or  0,  and  A  is  an  N-by-E  node-arc  incidence 
matrix  in  which  there  are  exactly  two  nonzero  entries  per  column, 
a  1  and  a  -1.  Any  convex  separable  quadratic  network  flow  pro¬ 
blem  with  N  nodes  and  E  edges  can  be  solved  by  such  an  LCP.  The 
major  reason  we  choose  this  LCP  in  our  numerical  study  is  that 
when  the  heuristic  approach  discussed  in  the  previous  section 
decides  to  set  some  x^  or  w^  to  zero,  it  is  very  easy  to  do  so 
to  reduce  the  problem  size.  For  more  general  problems,  it  may 
need  more  computational  effort  and  storage  space  in  doing  that. 
Beyond  that,  we  see  no  significant  difference  for  the  iterative 
method  and  its  heuristics  in  solving  LCPs  with  general  sparse 
positive  semidefinite  matrices.  So  the  following  results  should 
reflect  the  performance  of  the  method  not  only  for  the  special 
problem,  but  also  for  the  general  problems. 
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Many  problems  have  been  tested,  we  include  here  output  of 
five  runs  of  different  problem  size  ranges  and  different  methods 
for  choosing  the  tolerance  in  solving  (11).  The  following  remarks 
explain  what  each  column  means.  See  Table  1-5. 

Remark  5.1.  (Column  1)  All  problems  are  solved  in  less  than 
seven  iterations . 

Remark  5.2.  (Column  2)  As  stated  in  Section  2.1,  the  stopping 
criterion  is  that  |f(x)|  be  smaller  than  a  fixed  tolerance,  called 
Initial  Tolerance  in  each  output,  for  which  we  use  |1  q  1 1  /lOOOn 
and  find  it  satisfactory.  The  initial  tolerance  is  also  used  in 
solving  (11)  for  the  first  y.  In  all  except  Table  2  for  which 
the  same  tolerance  is  used  in  each  iteration  in  solving  (11), 
the  tolerance  is  halved  after  each  iteration.  One  extra  iteration 
is  done  to  bring  the  infeasibility  within  "user  specified" 
tolerance,  called  Final  Tolerance  in  the  outputs. 

Remark  5.3.  (Column  3)  The  program  is  given  the  solution  so  that 
it  can  compute  the  square  of  the  relative  error  ||  x^-x*||2/  H  x*||2 
We  choose  to  avoid  computing  the  square  root.  Note  that  the 
solution  procedure  is  in  no  way  dependent  on  the  given  solution. 

In  fact ,  the  time  for  computing  this  column  should  not  have  been 
counted. 

Remark  5.4.  (Column  4)  In  solving  (11),  a  cycle  is  complete  when 

all  of  the  n  constraints  w.>0  and  the  cut  have  been  checked  and 

1 

relaxed  if  violation  is  found. 
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Remark  5.5.  (Colvunn  5)  This  is  the  total  number  of  changes  of  y 
in  solving  (11).  Therefore  it  is  also  the  total  nximber  of  times 
of  all  constraints  violated. 

Remark  5.6.  (Column  6)  This  is  the  clock  that  starts  at  00:00.00 
in  each  run.  The  first  two  digits  are  in  minutes,  the  last  two 
are  in  hundreds  of  seconds . 

Remark  5.7.  Table  2  tells  us  that  high  accuracy  cannot  be  achieved 
without  reducing  the  tolerance,  which  is  consistent  with  theoreti¬ 
cal  results  of  (Ref.  14).  But  when  the  tolerance  is  reduced 
iteratively  in  Table  3 ,  the  computing  time  is  increased  greatly,  so 
is  the  number  in  Column  5.  The  latter  means  that  many  constraints 
are  violated  many  times ,  but  then  many  w^  and  will  be  set  to 
zero  by  the  heuristic  enhancement  criteria.  The  problem  size  can 
therefore  be  reduced  rapidly,  that  is  reflected  in  Table  4  in  which 
the  computing  time  is  less  than  one  percent  of  that  in  Table  3. 

The  same  remarkable  performance  is  found  in  Table  5 ,  in  which  a 
problem  of  the  total  of  116  variables  is  solved  with  relative 
accuracy  of  10  in  just  32  seconds--by  a  personal  computer. 

6 .  Additional  Observations  and  Summary 

In  all  our  test  problems  of  size  more  than  50,  starting  at  the 
second  iteration  all  constraints  violated  more  than  5  times  turned 
out  to  be  active,  and  more  than  90%  of  the  active  constraints  are 
violated  more  than  50  times  in  a  single  iteration.  In  other 
words,  the  program  has  picked,  without  a  single  error,  more  than 
90%  of  the  active  constraints  in  just  three  iterations.  Assuming 
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the  same  performance,  the  method  would  solve  problems  of  millions 
of  variables  in  minutes  or  seconds  on  mainframes  or  supercomputers. 
The  related  problem  of  how  to  reduce  the  problem  size  for  general 
LCP  while  keeping  as  much  sparsity  as  possible,  when  many  of  the 
active  constraints  are  known,  becomes  a  very  interesting  problem 
which  we  leave  to  the  future. 
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Table  1.  Output  By  LCPSOR,  no  heuristic  enhancement. 


Square  of  Number  of  No.  of  SOR  Clock 


Iteration 

f  (X) 

rel.  error 

cycles 

iterations 

0 

67.136070 

2.808698E-01 

3 

13 

00:00.05 

1 

10.386370 

i.373864E-02 

5 

23 

00:00.16 

2 

0.608786 

8.386110E-05 

17 

112 

00:00.49 

3 

-0.000018 

2.790467E-09 

28 

190 

00:01.04 

4 

-0.000044 

1.886973E-13 

53 

321 

00:02.09 

Number  of  nodes  =  3 ;  Number  of  edges  =  4 . 

Initial  tolerance  =  0.004000;  Final  tolerance  =  0.000005  . 
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Table  2.  Fixed  tolerance,  no  heuristic  enhancement. 


Square  of  Number  of  No.  of  SOR  Clock 


Iteration 

f  (X) 

rel.  error  cycles 

iterations 

0 

2369.166000 

1.815431E-01 

4 

46 

00:00.33 

1 

348.715500 

1.577661E-02 

6 

92 

00:01.10 

2 

19.965990 

2.935577E-03 

63 

848 

00:07.36 

3 

3.917515 

7.512022E-04 

314 

4235 

00:37.79 

4 

0.291592 

1.806640E-04 

747 

7048 

01:46.11 

5 

0.108357 

1.318387E-04 

280 

1678 

02:10.45 

6 

-0.082139 

1.292340E-04 

25 

129 

02:12.75 

7 

-0.082139 

1.292340E-04 

0 

0 

02:12.92 

Number  of 

nodes  =  20; 

Number  of  edges 

=  38. 

Fixed  tolerance  =  0.012897 
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Table  3.  Reducing  tolerance,  no  heuristic  enhancement. 


Iteration 

f  (X) 

Square  of  Number  of 

rel.  error  cycles 

No.  of  SOR 
iterations 

Clock 

0 

2369.166000 

1.815431E-01 

4 

46 

00:00.33 

1 

349.100200 

1.576520E-02 

6 

94 

00:01.10 

2 

20.196930 

2.908019E-03 

.  74 

1211 

00:08.73 

3 

4.116029 

7.416783E-04 

513'  • 

8657 

01:00.91 

4 

0.738284 

7.504172E-05 

1807 

32150 

04:06.95 

5 

0.106307 

2.941542E-06 

4287 

69869 

11:19.48 

6 

0.004395 

8.453883E-08 

5070 

70660 

19:25.02 

7 

0.000147 

7.029873E-11 

10185 

195481 

37:23.54 

Same  problem  as  in  Table  2 . 

Initial  tolerance  =  0.012897;  Final  tolerance  =  0.000005 
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Table  4.  Output  By  LCPSOR  with  heuristic  enhancement. 


Square  of  Number  of  No.  of  SOR  Clock 


Iteration 

f  (X) 

rel.  error 

cycles 

iterations 

0 

2369.166000 

1.815431E-01 

4 

46 

00:00.44 

1 

349.100200 

1.576520E-02 

6 

94 

00:01.43 

2 

20.196930 

2.908019E-03 

74 

1211 

00:10.77 

3 

3.963225 

4.940340E-04 

18 

82 

00:11.26 

4 

0.338470 

2.348035E-06 

39 

216 

00:12.25 

5 

-0.002658 

8.072329E-11 

93 

434 

00:13.73 

6 

-0.000023 

1.349873E-14 

98 

422 

00:14.21 

Same  problem,  same  initial  and  final  tolerance  as  in  Table  3. 
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Table  5.  Output  By  LCPSOR  with  heuristic  enhancement. 


Square  of  Number  of  No.  of  SOR  Clock 
Iteration  f(x)  rel.  error  cycles  iterations 


0 

3990.690000 

1.634287E-01 

8 

113 

00:01.32 

1 

627.943400 

1.455267E-02 

6 

159 

00:03.19 

2 

40.228890 

3.347591E-03 

71 

2540 

00:21.31 

3 

7.471679 

2.231230E-04 

14 

130 

00:21.97 

4 

0.767882 

2.733713E-06 

36 

436 

00:23.51 

5 

0.010428 

5.249066E-09 

123 

1146 

00:26.75 

6 

0.000015 

4.992850E-14 

200 

1784 

00:31.91 

Number 

of  nodes  =  40; 

Number  of  edges 

=  76. 

Initial 

tolerance  =  0 

.011388;  Final 

tolerance 

=  0.000005 
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